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AN EARTH ALBEDO MODEL 


Introduction 

Coarse Sun Sensors are often used on spacecraft as part of. Attitude 
Determination Systems. They function essentially as "Direction Cosine sensors, 
with their output approximately proportional to the cosine of the angle between 
the "boresight" of the sensor and a vector from the spacecraft to the sun. Their 
field of view is approximately hemispherical. 

The sensing element is typically a silicon solar cell which produces a current 
proportional to the energy flux incident on its surface. Energy can arrive however 
from sources other than the sun which lie within the field of view of the sensor. 

For an Earth-orbiting spacecraft, additional currents could be produced by such 
things as the Moon, sunlight reflected from some part of the spacecraft and/or 
radiation from portions of the Earth's surface. 

This report presents a mathematical model for the Coarse Sun Sensor output 
due to radiation originating from the sunlit portion of the Earth within the sensor 
field of view. 

Albedo 


Solar radiation arriving at the Earth's surface is generally considered to be 
partially absorbed, partially specularly reflected and partially diffusely reflected. 
Local surface characterisics and cloud cover conditions determine the relative 
importance of these phenomena. 

The energy which is absorbed is eventually re-radiated into space at infrared 
wavelengths. Solar cells are insensitive to this radiation. 

With specular reflection (as commonly occurs with mirrored surfaces) some or all 
of the incoming solar rays are reflected with the angle of reflection equal to the 
angle of incidence. Since a spacecraft would receive very little energy from even 
an entire Earth which was specularly reflecting this type of reflection is ignored 
here. 

Here, we consider the sunlit potion of the Earth to be a uniform, diffuse reflector 
and will use the word "albedo" in a limited sense, i.e. the albedo constant will be 
taken to be the ratio of the energy diffusely radiated from a surface to the total 
energy incident on the surface. 
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Diffuse reflectance 


Diffuse reflectance is due 
to the scattering of the 
incident light in all 
directions. Consider a 
small area of the Earth's 
surface, dAe- As the 
diffuse reflectance 
scatters, it forms a 2n 
steradian solid angle with 
constant intensity levels 
as shown. As the light 
travels outward, it spreads 
out, reducing the intensity. 


Incident Light 



2 n steradian solid 
angle of diffuse reflection 


ant intensity levels 


dA e , piece of sunlit earth in 


To determine the strength 

of the diffuse reflectance at the spacecraft's position, consider the 
geometry shown; 



where, n e : unit vector normal to dAe 

s : unit vector from the 
Earth to the sun 
B : point where s intersects 
the Earth's surface 


According to Wertz, F S un. the solar constant in the vicinity of the Earth, is 
approximately 1358 wat % 2 - The sunlight strikes the Earth with this 
intensity at point B. At locations away from this point, the intensity of the 
incoming sunlight decreases proportional to cosy, so that the solar flux 
reaching any given incremental area is: 

F* = F^„(n.-s) *•%, 

This incoming solar flux is partially absorbed and partially reflected. The 
amount of light reflected is proportional to the incident light by an albedo 
constant, ALB, which depends on the Earth's surface characteristics. (See 
Appendix II.) This model assumes that the albedo constant does 
not vary over the Earth's surface, neglecting the variation of 
diffuse reflectance with geographical features. A good estimate of 
the Earth’s annual average albedo constant is 0.3. 

Under these assumptions, the amount of solar flux diffusely reflected from 
dAe is given by: 

F out = ALB(F in ) 

= ALB(F sun )(n e s) " at %, 2 


2 


As the light travels outward, it spreads out, reducing its flux, so only a 
portion of the diffuse reflectance from dAe actually reaches the 
spacecraft's position. Consider a hemisphere centered at dA© with radius 
D, the distance from dAe to the spacecraft's position. 


/ 1 

/ !© - 


'yj o) 0 : maximum flux at distance D 
r- ^ Spacecraft 

A-v^. flux at spacecraft: 


CO a |b = t0 o C0S© watts/m 2 


■—I i.Jl 


dAe 


total energy over this hemisphere = J co 0 cos<|>(2jtDsin<t>)Dd<t> 


0=0 


= 7cD 2 co 0 watts 


By conservation of energy, the total energy over this hemisphere must 
equal the energy radiated from dAe. Using this fact, solve for ©o : 


tcD 2 co 0 = F 0Ut dA e 

_ A lB(F sun )(n e s)dA e wa tts/ 

0)0 " tiD 2 /m 2 

The solid angle subtended by dA e from the spacecraft's point of view, is 



So, the flux at the spacecraft's position is: 


©alb = ©o cosQ 

ALB(F sun )(n e s)dA 

©alb “ Z 


watts/ 

/m 2 
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The diffuse reflectance from dAe only affects the output current of the sun sensors 
if the light reaches the sun sensors. For this to occur, the area, dAe must meet 
the following conditions: 


D 


Assume a perfectly spherical Earth. The area, dAe, must be on the 
sunlit side: |(s n e ) >~0 



The area must be in the spacecraft's field of view: (a < a max ) 



Spacecraft 


3) Assume that a Coarse Sun Sensor (CSS) has a conical field of 
view with a half angle of A. The reflected light must be in the CSS 
field of view: (u • ness) - cos A 
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Sun Sensors 

This algorithm models Coarse Sun Sensors. In general, the Coarse Sun 
Sensor's output current, less, is proportional to cosp, where p is the angle 
between the CSS boresight and u, the direction of light source: 

n egg light source less 06 (^css^OSP 

where, co^ : input flux to CSS 

n css : unit vector normal to the CSS 

u : unit vector from the CSS to the 
css light source 



Coarse sun sensors typically produce output currents (in the microampere range) 
when illuminated by direct and/or reflected solar energy. In this 
(!) analysis, the output current (less) is normalized such that its value 
Tincoming is 1-0 when the sun ,ies alon 9 the boresight of the CSS and there 
y sunlight is no additional heat input from any other source. 

"css A 

■1 M0 

u'P 

CSS 


In reality, the CSS receives incoming light from Earth as well as from the sun. 

Consider a small area of the Earth as the sole light source. 
\dA e in this case, the incoming light is albedo rather than 

"css VAfoedoSunlight, so cD css =cD a ib, the energy flux due to albedo at the 
A r spacecraft's position. 



css 


If sunlight along the boresight produces a current of Imax. the maximum CSS 
current due to albedo is: 


imax _ 
^sun 

dlalbmax 


dlalbmax 

©alb 

l max (ALB)(n e s)dA 

n 


5 


Since, less is proportional to (5, use the maximum values to find an appropriate 
proportionality constant for dlalb. the CSS output current due to albedo. 

d^alb = ^( w alb )COSP 

Substitute the maximum values to find K. Assume that the direction of the 
incoming light coincides with the CSS boresight such that (5 = 0: 

(l max )(ALB)(n e s)dA = j ALB(F 8un )(h e s)dA ~ 
n |_ n 

K = U. amperes 

p waits/ 

'sun /m 2 

In general, CSS output current is calculated: 

Us = K(co)cosp 

Note that cosp = (n css u); substitution into the general equation, gives theCSS 
current due to albedo from dAe as: 
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Calculations 

Consider the Earth as a light source. To calculate dl a ib. it is necessary to know u, 
the direction of the incoming light. At great distances, the Earth may be 
considered a point source; however, at orbital altitudes the curvature of the Eart 
must be taken into account. To accomplish this, the Earth is divided into 
incremental areas; each reflects light at a slightly different angle. The CSS 
current due to albedo also depends on the strength of the diffuse refiectance 
which is a function of ALB and Fj n . As previously discussed, it is necessary to 
define n e , a unit vector normal to a given area if the Earth, to calculate the 
strength of the incoming solar flux. 



Earth 


The incremental areas are found by dividing the solid angle subtended by t 
Earth from the spacecraft's point of view into NA solid angles of equal size (dA), 
each one corresponding to a piece of the Earth (dAe). Once the areas are 
defined, the necessary unit vectors are calculated for each area. 


Dividing th e Field of View 

If the spacecraft is in a circular orbit, the size of the spacecraft s 
field of view remains constant, so that once these areas are defined they 
are fixed. If a m ax is the angle from the subsatellite point to the horizon for 
a spacecraft at altitude. ALTP, above a spherical earth of radius R e , then: 



From this, the solid angle subtended by the earth is calculated: 
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A max = 2lt(1 ’ C0SCt max)- 


Consider a sDhere of unit radius centered on the spacecraft, where the 
soTd angte subtended the Earth defines a spherical segment with a 

surface area of Amax: 



Earth 


To divide the field of view into NA equal parts, a set of N concentric -circles 
(i= 1 , N ) is drawn on the spherical segment, where the largest circle 
eorresDonds to the horizon of the Earth. Each circle defines anot 
spherical sub-segment and represents the intersection of a cone with the 

unit sphere. 

Next, each sub-segment must be subdivided into equal solid angles. 

To understand how this system of concentric circles may be subdivided 
into NA equal areas, consider a flat circle: 



8 









Define the centers of these areas as; one at the center °^ e ^ c '® 

and eight others spaced around the dashed circle which divides eac 
Ire a into two equal halves; the centers are located at the points where 
this dashed circle intersects the radial lines which bisect each area. In this 
case, the dashed circle must enclose an area of 5 units. Let r 2 be the 
radius of the dashed circle. Thus, 

jrr 2 = 57tr? 

r 2 = r-,V5 

Thus, its radius will be V5 times that of the original unit circle. 

In general, let : N = total number of circles 

NA = total number of equal areas 

Rj = relative radius of i th circle 

rj = relative radius of circle locating the 

centers of the i th annulus areas 

0j = angle between radial dividing lines of the i^ circle 

Next add a circle of relative radius Rj = 5. The total area of the new 
annulus formed will be 25 - 9 = 16 units. Divide this annulus into 16 equal 
parts. Now, the figure will have 25 equal areas. The centers of these new 

areas will lie on a circle of radius, rj =V 17 

Continuing to add circles which form equal width annuluses and dividing 
each into eight more sections than the previous circle, the result is. 


In general, 


i 

Bi 

ii 


2 

3 

V5 

45.0° 

3 

5 

Vl7 

22.5° 

4 

7 

V37 

15.0° 

5 

9 

V 65 

11.25 

6 

11 

VIM 

9.0° 

NA = 

(2N-1) 2 



= 

^4(i-1) 2 +1 

, i = 

2 to N 


G. = 


360 

8(1-1) 


i = 2 to N 
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Now, extend this planar analogy to the spherical segment. Begin with the 

spherical segment with a half cone angle of 

ctmax, corresponding to the spacecraft's field 
of view. The surface area of this unit cone is 
equivalent to the solid angle subtended by 
the Earth, Amax- 

If this area is to be divided into NA equal 
parts, then the area of one unit is: 



A, = 


'max 

NA 


Let ai define a central sub-segment of area, 
A-j. Then, 

2jc(1 - cosa-i) = ^max 
NA 


cosa-i = 


NA - 1 + co soci 
NA 


max 


Define the corresponding unit vector, u 1t through the center of the area. In 

this case, u-| is along the axis of the field 
of view cone. 

Now find a cone which encloses an area 
of 5Ai; this is analogous to the first 
dashed circle, corresponding to i = 2: 



2tc( 1 - cosa 2 ) = 
1 - cosa 2 = 
cosa 2 = 


5 A, 


max 


NA 


— (1 - cosa max ) 
NA 

NA - 5 + 5cosoc 


max 


NA 




In addition to u-|, define eight more unit vectors 
at an angle of 012 from the axis of the field of 
view cone, equally spaced around the axis in 

increments of ©2. where ©2 = 

These are the first nine unit vectors. 


Next, find a cone which encloses an area of 17Ai; following the planar 
example, the centers of these areas will be evenly spaced around this 

circle in increments of ©3, where ©3 = l n this case, 


cosa 3 = 


NA-17 + 17cosa max 
NA 


resulting in 16 unit vectors at an angle a 3 from the axis. 


In general, there are NA solid angles of equal size. To define the 
appropriate unit vectors it is necessary to calculate the location of the 
center of each dA. 


For the k th solid angle (k=1,NA), the center's coordinates are (Yk,5k). 
where y k is the angle from the axis of the unit field of view cone to the 
center of the k* h dA and 5 k is the radial angle measured counterclockwise 
from the 1 axis to the center of the k th dA. 

The angular coordinate is calculated as follows: 

Tj = distance to center of k th area 
(constant for each circle) 

= -^4(i - 1)^ + 1 



yk coordinate of k*h area is: 


cosYk =cosoij = 


NA - r? + tf (cos a ma x) 


NA 


for 


all areas in the i tn circle, where i = 2 to N 




To calculate the size of the radial division of the i th circle, let; 



circle. 


number of radial divisions in 

i th circle 
8(i - 1) 

size of radial division in i th 
circle 

360° 

M 


where j = 1 to M for the i^ 


Orbital Coordinate Frame 

The unit vectors are calculated in a 1-2-3 orbital coordinate system, which 
moves with the spacecraft. If the orbit is circular, the vectors are fixed 
within this frame, and n e and u need only be calculated once, then 
transformed to the spacecraft body frame at each time step. 

The axes of this 1-2-3 frame are defined as: 



where, R : spacecraft unit position vector in earth centered coordinates 
V : spacecraft velocity unit vector 
H 0 : unit orbit normal vector 
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Calculating u 

In the 1-2-3 coordinate frame, a unit vector from the spacecraft to dAe is: 


Uk 


siny^ cos 5k 
siny k sin5 k 


-cosYk 


J123 frame 


Calculating 

Once u has been found for each area, calculate n e , the normal to each 


area. First, consider the 1-2-3 coordinate frame. In general, any point G in 
this coordinate system may be located by G 123 , a vector from the Earth's 
center to the point, G. 


A 

1 



a point on the Earth’s surface 


A 



Gl23 

g 123 


Du, 

Du 2 

R + Du 3 
R + Du 


where D is the distance from the 
spacecraft to the point G. 

To define the normal vector, G must be 



^/(Dui ) 2 + (Du 2 ) 2 + (Du 3 + R) 2 = R e 

Solve this quadratic for D: 

D = -Ru 3 ± V(Ru 3 ) 2 - (R 2 - R|) 


Choose the smaller value of D, corresponding to G, not G'. 
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Once the appropriate D is calculated, the normal unit vector is: 



Implementing the A lgorithm 

Appendix I contains a copy of the FORTRAN subroutine based on this paper. 
The table of variables provides a cross reference between the notation in the 
paper and the subroutine variables. 


In order to optimize the calculations, the subroutine uses the conditions to 
determine which areas actually affect a given CSS, then calculates the output 
current based on only those areas. 


During the first call, the subroutine initializes constants and calculates the unit 
vectors, u^and n ek , for each of the NA areas (k=1,NA). To calculate these 


vectors, begin with the first circle (i=1), and place u-j along the axis of the 
spacecraft's field of view cone. This corresponds to the vector (0,0,-1) in the 

1-2*3 orbital frame. 



For the second circle (i=2), calculate 0 C 2 
and 02 , the size of the angular and radial 
divisions for this circle. Then, find the 
coordinates (a^.S^) of the centers of the 
M areas comprising this circle. For each 
set of coordinates, calculate u^, in the 
1-2-3 orbital frame. Continue this 
process, until u^has been defined for 
each radial division in the second circle, 
(This corresponds to areas k=2,9). Move 
to the next circle and repeat this process 
until Ok has been found for all areas, 1 
through NA. Once all of the 0^ vectors 
are calculated, the corresponding 
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n e vectors are found. By defining the areas based on the spacecraft s field of 

view not only is Condition 2 automatically satisfied, but the smallest possible 
area of the Earth's surface is considered. This ends the initialization process. 

At each call, the subroutine calculates the sun unit vector in the 1-2-3 orbital 
frame; the position of the sun vector is necessary to test Condition 2. For the I 
CSS (1=1, P), use the remaining conditions to determine if the albedo from the k 
area effects this CSS. First, determine if the k* h area is lit (Conditionl ).lfso 
determine if it falls within this CSS's field of view (Condition 3). Ne>d, check that 
the incoming light is not blocked. (Refer to the Operational Notes in Appendix \ for 
a complete explanation of blockage.) If all of these conditions are met, calculate 
dlalb from this area. Once dl a ib has been calculated for areas 1 through NA 
sum these currents to calculate less* the total CSS output current due to albe o 
for the |th sun sensor. Repeat this process for each CSS, 1 through P. 


Conclusions 

This simplified albedo model was developed for use in spacecraft control system 
simulations, specifically, for modeling Coarse Sun Sensors. It is based on 
several approximations. Only diffuse reflectance is included; specular 
reflectance is neglected. For an elliptical orbit, the unit vectors associated with 
the incremental areas should change direction with altitude; instead, this 
alaorithm assumes a circular orbit. The albedo constant is set to the annual 
qlobal average for the entire Earth; Appendix II illustrates how the percentage of 
light reflected truly varies with geographical features. The Earth is considered a 
perfect sphere which does not rotate; it was unnecessary to model rotation since 
the albedo constant was not varied. 
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Table of Variables 


Units are given in parentheses. 

The equivalent variables from the albedo subroutine are given in 

brackets. 


A : transformation matrix from spacecraft body coordinates to the 1-2-3 orbital 

frame [A] 

Amax : solid angle subtended by the portion of the Earth's surface, visible from the 

spacecraft (steradian) [AMAX] 

Ai: area of central subsegment of unit sphere, define as one unit area, equal to 

the size of dA (steradian) [A!] 

ALB : ratio of diffuse reflectance to incident light [ALB] 

ALTP : altitude of the spacecraft at perigee (km) [ALTP] 

B : point where s intersects the Earth's surface 

BLOCK : a Px4 matrix where the first three elements are the components of the 

unit vector along the axis of the blockage cone in spacecraft 
body coordinates, and the fourth element if the half angle of 
the blockage cone in degrees [BLOCK] 

BLOCKAGE : flag used to indicate whether the Coarse Sun Sensor’s field of view 

is partially blocked [BLOCKAGE] 

D : distance from the k^ dAe to the spacecraft's position (km) [D] 

dA : solid angle subtended by the area dA e in the spacecraft's field of view 

(steradian) 

dAe : a piece of the sunlit Earth's surface (m2) 

dlalb : normalized Coarse Sun Sensor output current due to the albedo from a 

single dAe (amperes) [Dl] 

dlalbmax : maximum possible Coarse Sun Sensor output current due to the 

albedo from a single dAe (amperes) [IMAX] 

Fjn : solar flux input to a given dAe ( wa, %,2) 

Fout : solar flux output from a given dAe ( wat %,2 ) 

Fsun : solar constant in the vicinity of the Earth ( wat %,2) 

G : any point in the 1-2-3 orbital coordinate frame 

G 123 : a vector from the Earth's center to the point, G, in the 1-2-3 orbital frame 

H 0 : unit orbit normal vector [HO] 
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•max : maximum Coarse Sun Sensor output current, (amperes) [CSSIMAX] 
less : Coarse Sun Sensor output current (amperes) [ALBICSS] 
j : refers to the jth radial area in the i*h annulus [J] 

K : proportionality constant for normalized Coarse Sun Sensor output current due 

to the Earth's albedo [KAPPA] 

k : refers to the k*h unit area out of NA equal areas [K] 

I : refers to the l*h Coarse Sun Sensor of out P Coarse Sun Sensors 

M : number of radial division in the i th annulus circle [M] 

N : number of concentric circles [N] 

NA : number of equal areas Earth is divided into [NA] 

n css : unit vector normal to the Coarse Sun Sensor [NCSS] 

n e : outward pointing unit vector normal to dAe [NORM] 

P : number of Coarse Sun Sensors [P] 

Rj : relative radius of i*h circle 

R : unit vector from the Earth's center to the spacecraft's position [RB] 

R : vector from the Earth's center to the spacecraft's position (km) [RMAG = |r|] 
Re : radius of the Earth (Km) [RE] 

H : perpendicular distance from the axis of the spacecraft field of view cone to the 

center of the i th annulus; given in terms of ri [Rl] 

ri : radius of central circle of unit area 

s ; unit vector from the Earth to the sun [SB] 

0 : in general, a unit vector from the Coarse Sun Sensor to the source of the 

incoming light; specifically, it is a unit vector from the 
spacecraft to the k th dAe [U] 

V : spacecraft velocity unit vector [VB] 

otj : angle from the axis of the unit radius spacecraft field of view cone to the center 

of the i*h annulus 

otmax : half angle of the cone encompassing the spherical segment; the angle 

from the nadir to the horizon [ALPHAMAX] 

ai : half angle of cone encompassing a central sub-segment of unit area 

3 : angle between u and n css 

A : half angle of Coarse Sun Sensor's conical field of view [CSSLM = cosA] 
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5k : radial coordinate of the center of the k th dA; the angle measured 

counterclockwise from the 1 axis to the center of the k* dA 

7k • angular coordinate of the center of the k th dA; the angle from the axis of the 

unit radius spacecraft field of view cone to the center of k™ 
area; it is constant for all areas in the ft* circle [CGAMMA 
and SGAMMA, cosine and sine of 7k] 

© ; angle between the k^ h e and the spacecraft's position 

0j : size of the radial division of the i th circle [THETA] 

<j> : variable of integration 

©alb : flux at the spacecraft’s position ( wa % 2 ) 

©css ; flux input to the CSS ( wat % 2 ) 

©o : maximum flux at distance D ( wat %, 2 ) 
y : angle between s and the k th h e 
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Operational Notes 

Before running the subroutine the user MUST set the following constants: 
BLOCKAGE : 0 = no blockage; 1 = blockage exists 

N : number of circles. DO NOT EXCEED N=16 or all arrays currently 
dimensioned to 1000 must be redimensioned to NA, where NA - 

(2N- I) 2 

CSSLM : cosine of the half angle of the Coarse Sun Sensor field of view 
(cosA) 

P : number of CSS. DO NOT EXCEED 8 or all arrays currently 
dimensioned to 8 must be redimensioned to P. 

ALTP : Altitude of the spacecraft at orbit perigee in kilometers. 


ALB : albedo constant 

CSSIMAX : CSS current with full Sun directly along CSS boresight 

All of these constants, except BLOCKAGE, are used are previously defined 
BLOCKAGE is a switch, indicating whether the CSS's field of view is partially 
blocked. This subroutine models blockage as a cone, as illustrated below. 

If blockage exists, the locations must be listed in the BLOCK matrix wN<J is 
passed in at each call. If the CSS are not blocked, simply set the BLOCKAGE 
flag to zero, and zero the BLOCK matrix. 

The BLOCK matrix is a Px4 array where P is the number of CSS. 



The first, second, and third elements are 
the components of a unit vector along the 
axis of the blockage cone in spacecraft 
body coordinates. The fourth element is 
X, the half-angle of the blockage cone: 

CSS Xs/c Ys/c Zs/c ^deg 

1 

2 


P 


The following variables are passed in at each call: 

T : time 

NCSS : CSS boresight unit vectors in spacecraft body coordinates 
BLOCK : array of CSS blockage values 

VB : spacecraft velocity unit vector in spacecraft body coordinates 
RB : spacecraft position unit vector in spacecraft body coordinates 
SB : sun unit vector in spacecraft body coordinates 


PAGE < INTENTIONALLY BLANK 2 5 


PPECSDtNS P5Q? PLANK NOT Him 


The albedo subroutine returns 

ALBICSS : final CSS currents due to albedo (in same units as CSSIMAX) 


26 


Albedo Subroutine 


SUBROUTINE ALBEDO(T,NCSS,ALBICSS, BLOCK, VB.RB, SB) 


* This subroutine approximates the current in Coarse Sun Sensors due 

* to the Earth's albedo. Only includes albedo due to diffuse reflectance. 

*** ASSUME: 

* neglect albedo due to specular reflectance 

* circular orbit 

* perfectly spherical earth 

* conical field of view for coarse sun sensors 

*** CONSTANTS: 

* ALPHAMAX : S/C field of view half angle 

* AMAX : solid angle subtended by the portion of the Earth's 

* surface visible from the spacecraft (steradian) 

* DA : solid angle subtended by dAe, a small piece of sunlit 

* earth, from the spacecraft's field of view (steradian) 

* RE : radius of the earth (km) 

* RMAG : distance from center of earth to S/C in circular orbit 
** User sets these constants in the ALBEDO subroutine: 

* ALB : albedo constant: ratio of diffuse reflectance to incident 

* light 

* ALTP : altitude of spacecraft at perigee (km) 

* BLOCKAGE: flag to indicate whether the CSS are blocked. [0.0 = no, 

* 1.0 = yes] If they are the blockage numbers must be in 

* the (P,4) array BLOCK 

* CSSIMAX : maximum current output by CSS to Attitude Control System 

* when sun lies along CSS boresight. (Include any scaling by 

* Attitude Control Electronics.) 

* CSSLM : cosine of the half-angle of the CSS FOV; ALBEDO assumes 

* all of the P CSS have the same size conical FOV 

* N : number of concentric circles 

* NA : number of equal areas 

* P : number of coarse sun sensors 

** User passes in these constants from the calling routine at each time 
** step: 

* NCSS : unit vectors of CSS boresight in S/C body coordinates 
*** VARIABLES: 

** User passes in these variables at each time step: 

* BLOCK : a Px4 array which holds the blockage values for the CSS 

* Modeled as a blockage 'cone'. First three elements of each 

* row are the S/C body coordinates of this cone's axis. 

* The fourth elements is the angle of the blockage cone 

* in degrees. 

* RB : spacecraft position unit vector from earth center to S/C in 

* S/C body coordinates 

* SB : unit vector from the Earth to the sun in S/C body coordinates 
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T : time 

VB : spacecraft velocity unit vector in S/C body coordinates 
Program variables: 

A : 3x3 transformation matrix from S/C body coordinates, to 1-2-3 frame 
A1 : area of central subsegment of unit S/C centered sphere, 
defined as one unit area equal to the size of DA 
BLOCKED : test flag. If the Lth CSS blockage blots out the 
albedo light from the Kth infinitesimal area, 

BLOCKED = 0.0 (zeros the Dl from that area). If not 
blocked, BLOCKED = 1.0 

CGAMMA and SGAMMA : cosine and sine, respectively, of the angle 

from axis of FOV (Field of View) cone to normal 
of Kth dAe (deg) 

CSSIMAX : maximum possible CSS output 
D : distance from Kth dAe to the spacecraft's position 
Dl : CSS output current due to albedo from Kth dAe 
HALFTHETA : equals THETA/2 

HO : unit vector in direction of orbit normal; unit vector for axis 2 
I : refers to lth concentric circle out of N concentric circles 
J : refers to Jth radial area in lth annulus 
K : refers to Kth unit area out of NA equal areas 
KAPPA : reflective constant for Kth dA 
L : refers to the Lth CSS out of P CSS 
LIT : flag, indicates if Kth dA is in light or darkness 
M : number of radial divisions in lth circle 
NDOTU : CSS head unit vector dotted with Kth S/C to ground unit 
vector to determine if light from Kth dAe is visible 
to the spacecraft 

NORM : outward pointing normal vector for Kth dAe 
ONE : HO x RB; unit vector along axis 1 

Rl : distance from the spacecraft FOV cone's axis to area centers for lth 
annulus 

SI 23 : sun unit vector in 1-2-3 frame coordinates. 

THETA : size of radial divisions in lth circle (deg) 

U : unit vector from S/C to 'center' of Kth dA© (location of Kth 
dA's NORM) in 1-2-3 coordinates 
U123 : temporary holding vector for each U 

USC : U transformed from 1-2-3 coordinates frame to S/C body coord. 
OUTPUT: 

ALBICSS : final CSS currents due to albedo 


IMPLICIT REAL*8 (A-H.O-Z) 

INTEGER l,J,K,L,M,N,NA,NATH,P,KL 

REAL*8 T, RE, ALTP.ALPHAMAX, PI, RMAG, HALFTHETA, AMAX, CSSIMAX 
REAL*8 NDOTU, CSSLM, THETA, Rl, CGAMMA, SGAMMA, ALB.A1 
REAL*8 UDOTBLOCK, BLOCKAGE, BLOCKED 
REAL*8 Dl, KAPPA, SDOTN,SB(3), Si 23(3), A(3,3),VB(3),RB(3),HO(3) 



REAL*8 ONE(3),U123(3) I USC(3),BLOCK(8, 
REAL*8 U(1000,3),NORM(1000,3),D(1000) 
REAL*8 LIT(1000),INALBED 
COMMON/ALBFLAG/INALBED 


4) 

,NCSS(3,10),ALBICSS(1 0) 


begin initialization loop: 

IF (T .EQ. 0.0) THEN 
set constants: 

BLOCKAGE = 1.0 
N = 6 

NA = (2*N-1)**2 
P = 8 

RE = 6378.0D0 

CSSLM = DCOSD(SO.ODO) 

ALTP = 350.0 
RMAG = RE + ALTP 
PI = DACOS(-I.ODO) 

ALPHAMAX = DASIND(RE/RMAG) 

AMAX = 2*PI*(1 -DCOSD(ALPH AMAX)) 

A1 = AMAX/NA 

ALB = 0.30D0 

CSSIMAX = 1.0D0 

KAPPA = ALB*CSSIMAX*A1/PI 

calculate unit vectors from spacecraft to Kth area normal 
K = 1 
DO 1=1, N 
M = 8*0-1) 

IF (M .EQ. 0.0) THEN 
U(K,1) = 0.0D0 
U(K,2) = 0.0D0 
U(K,3) = -1.0D0 
K= K+1 
ELSE 

THETA = 360.0D0/DBLE(M) 

HALFTHETA = 0.5D0*THETA 
Rl = SQRT(DBLE(4*(I-1)*(I-1)+1)) 

CGAMMA = (NA - Rl**2 + RI**2*(DCOSD(ALPHAMAX)))/NA 
SGAMMA = SQRT(1 -CGAMMA**2) 

DO J=1,M 

U(K,1) = SGAMMA*DCOSD(THETA*J - HALFTHETA) 
U(K,2) = SGAMMA*DSIND(THETA*J - HALFTHETA) 
U(K,3) = -CGAMMA 
K= K+1 
END DO 
ENDIF 
END DO 
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* calculate area normals 
C = RMAG**2 - RE**2 
DO K= 1,NA 

B = RMAG*U(K,3) 

D(K) = -B-SQRT(B**2-C) 

NORM(K,1) = D(K)*U(K,1)/RE 
NORM(K,2) = D(K)*U(K,2)/RE 
NORM(K,3) = (RMAG+D(K)*U(K,3))/RE 
ENDDO 

ENDIF 

* end initialization loop 

* AT ALL TIMES: 

* FIND matrix A to go from s/c body to orbital frame 1-2-3 
CALL CROSSU(HO,RB,VB) 

CALL CROSSU(ONE,HO,RB) 

DO KL=1,3 
A(1,KL) = ONE(KL) 

A(2,KL) = HO(KL) 

A(3,KL) = RB(KL) 

ENDDO 

* transform sun unit vector from s/c body to 1-2-3 
CALL T0123(S123,A,SB) 

* Find current in Lth CSS due to albedo of the lit areas: 

DO L = 1,P 

ALBICSS(L) = 0.0D0 
DO K=1,NA 

c determine if Kth area is lit: 

SDOTN = 

$ (S123(1)*NORM(K,1)+Sl23(2)*NORM(K,2)+Sl23(3)*NORM(K,3)) 

IF (SDOTN .GT. 0.0D0) THEN 
c determine is the Kth dA is the Lth CSS's FOV: 

DO KL=1,3 
U123(KL) = U(K,KL) 

ENDDO 

CALL TOSC(USC,A,U123) 

NDOTU=(NCSS(1 ,L)*USC(1)+NCSS(2,L)*USC(2)+NCSS(3,L)*USC(3)) 
IF (NDOTU .GT. CSSLM) THEN 
If there is CSS blockage, determine if it blocks 
the albedo from Kth area: 

BLOCKED = 1.0D0 
IF (BLOCKAGE .EQ. 1.0) THEN 
UDOTBLOCK=BLOCK(L,1)*USC(1)+BLOCK(L,2)*USC(2) 

$ +BLOCK(L,3)*USC(3) 
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IF (UDOTBLOCK .GT. DCOSD(BLOCK(L,4))) BLOCKED = 0.0D0 
ENDIF 

Dl = KAPPA*NDOTU*SDOTN*BLOCKED 
ALBICSS(L) = ALBICSS(L) + Dl 
ENDIF 
ENDIF 
ENDDO 
ENDDO 

RETURN 

END 


SUBROUTINE CROSSU(C,A,B) 

C = AxB; normalize C and return a unit vector 
INTEGER I 

REAL*8 C(3),A(3),B(3) 

C(1) = A(2)*B(3)-A(3)*B(2) 

C(2) = A(3)*B(1)-A(1)*B(3) 

C(3) = A(1)*B(2)-A(2)*B(1) 

CM AG = SQRT(C(1)**2+C(2)**2+C(3)**2) 

DO 1=1,3 

C(l) = C(I)/CMAG 

ENDDO 

RETURN 

END 

SUBROUTINE T0123(V123,A,VSC) 

transforms a vector from S/C body to 1-2-3 orbital frame 

INTEGER UK 

REAL*8 V123(3),A(3,3),VSC(3) 

DO J=1,3 
DO 1=1,3 

V123(J) = V123(J) + A(J,I)*VSC(I) 

ENDDO 

ENDDO 

VMAG = SORT (V 1 23(1 )**2+V 1 23(2)**2+V 1 23(3)**2) 

DO K = 1 ,3 

V123(K) = V123(K)/VMAG 
ENDDO 
RETURN 
END 

SUBROUTINE TOSC(VSC,A,V123) 

transforms a vectors from orbital frame 1-2-3 to S/C body 

INTEGER I.J.K 

REAL*8 V123(3),A(3,3),VSC(3) 

DO 1=1,3 
DO J=1 ,3 
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VSC(I) = VSC(I) + A(J,I)*V123(J) 

ENDDO 

ENDDO 

VMAG = SQRT(VSC(1)**2+VSC(2)**2+VSC(3)**2) 


DO K=1 ,3 

VSC(K) = VSC(K)A/MAG 
ENDDO 
RETURN 
END 
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